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COMBINED  SUMMARY 


Two  different  approaches  to  the  acceleration  of  iterative  algorithms 
for  the  numerical  solution  of  differential  systems  have  been  developed. 
General  form  of  the  non-linear  minimal  residual  method  has  been  analyt- 
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ically  determined  and  numerically  confirmed  foy  linear  and  non-linear 
problems.  The  method  was  applied  to  multi-step  algorithms  for  effectively 
determining  optimal  values  of  each  of  the  acceleration  parameters  at  each 
time  step.  It  was  found  that  both  the  rate  of  iterative  convergence  and 
the  smoothness  of  the  iterative  convergence  can  be  substantially  aug¬ 
mented  by  the  use  of  these  multiple  optimal  acceleration  parameters. 

The  second  approach  involves  a  composite  adaptive  method  which  is 
based  on  variational  techniques.  An  automatic  procedure  for  determining 
splitting  parameters  needed  in  the  iterative  solution  of  large  sparse 
linear  systems  was  developed.  It  was  then  complemented  with  the  gener¬ 
alized  conjugate  gradient  acceleration  procedures  and  successfully  ap¬ 
plied  in  the  symmetric  successive  overrelaxation  method  and  in  the  shifted 
incomplete  Cholesky  method. 

\  ”  1 

A.  Report  of  the  G rr.-.o  Headed  by  Dr.  G.S.  Dulikravich 

A .  1  Abstract 

A  generalized  non-linear  minimal  residual  (GNLMR)  method  based  on 
the  time-dependent  approach  and  multi-step  algorithm  has  been  developed 
to  accelerate  iterative  methods  for  solving  linear  and  non-linear  prob¬ 
lems.  Both  theoretical  studies  and  numerical  experiments  show  the 


monotone  convergence  behavior  of  the  GNLMR  method.  It  is  found  that  the 
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rate  of  convergence  and  the  smoothness  of  convergence  for  iterative 
methods  can  be  improved  by  the  optimized  multi-step  algorithm. 

A. 2  Research  Objectives  and  the  Statement  of  Work 

The  NLMR  method  developed  originally  by  S.  R.  Kennon  is  capable  of 
accelerating  Iterative  methods  for  the  solutions  of  linear  and  non-linear 
problems.  This  method  can  be  considered  as  a  time-dependent  technique 
with  variable  time  step  size.  The  basic  idea  of  this  method  Is  based  on 
the  sequential  minimization  of  the  global  residual.  A  relaxation  factor 
is  introduced  to  minimize  the  l_2  norm  of  the  residual  at  each  iteration. 
This  iteration-dependent  optimal  relaxation  factor  drives  the  iterative 
solution  to  convergence.  In  order  to  obtain  a  smooth  convergence  history, 
several  options  are  available  to  this  acceleration  procedure.  For  in¬ 
stance,  each  three  regular  iterations  can  be  followed  by  two  accelerated 
Iterations. 

Recently,  this  method  has  been  generalized  by  Huang,  Kennon  and 
Dulikravich  by  using  time-dependent  method  and  multi-step  algorithm. 
Both  theoretical  studies  and  numerical  experiments  proved  the  monotone 
convergence  behavior  of  this  method.  With  the  multi-step  algorithm,  it 
was  found  that  the  rate  of  convergence  and  the  smoothness  of  convergence 
of  the  NLMR  method  can  be  improved  even  further. 

Since  the  limitation  on  the  time  step  size  and  the  optimal  value  of 
the  time  step  size  can  analytically  be  determined  with  the  GNLMR  method. 
It  is  believed  that  Euler  and  Navier-Stokes  equations  of  gasdynamics  can 
be  efficiently  solved  using  the  GNLMR  method. 


A. 3  Status  of  the  Research 


A. 3.1  Introduction 

The  time-dependent  technique  proposed  by  Moretti  and  Abbett  [1]  in  the 
mid-1960's  was  the  first  successful  method  to  solve  problems  governed  by 
equations  of  mixed  type.  The  steady  state  solution  was  obtained  by 
starting  with  the  unsteady  equation,  and  marching  the  solution  along  the 
time  coordinate  until  convergence  was  achieved.  Nowadays,  the 
time-dependent  method  is  widely  used  in  computational  fluid  dynamics  for 
the  solution  of  the  steady  state  Euler  and  Navier-Stokes  equations  [2], 
[3],  [4]. 

It  is  interesting  to  note  that  most  iterative  methods  for  the  solution 
of  steady  state  problems  can  be  shown  to  be  equivalent  to  methods  for 
solving  time-dependent  problems  of  either  parabolic  or  hyperbolic  type 
[5],  [6].  The  relaxation  factor  used  in  accelerating  an  iterative  method 
to  obtain  the  converged  solution  plays  the  same  role  as  the  time  step  size 
in  advancing  the  transient  solution  to  the  steady  state  solution  for  a 
time-dependent  problem.  This  approach  offers  several  advantages. 

For  example,  the  mechanism  of  an  acceleration  scheme  can  be  unveiled 
and  an  optimal  value  of  relaxation  factor  (optimal  time  step  size)  could 
be  analytically  determined.  If  accurate  time  evolution  is  required  for 
an  unsteady  problem,  the  time  step  size  should  be  small  to  guarantee  both 
the  stability  and  the  accuracy  of  the  solution.  Consequently,  very  often 
a  compromise  must  be  made  beetween  the  computer  time  requirements  and  the 
accuracy  of  the  solution  of  such  computation  provided  the  method  is  sta¬ 
ble.  Since  the  limitation  on  the  time  step  size  can  be  analytically  de¬ 
termined  using  the  GNIMR  method,  a  satisfactory  compromise  can  be 
achieved.  If  transient  behavior  is  of  no  interest  to  us,  it  is  convenient 
to  interpret  the  physical  time  as  the  artificial  time  or  to  reformulate 
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the  governing  equations  in  terms  of  certain  artificial  time.  Then  the 
GNLMR  method  can  be  applied  to  determine  the  optimal  value  of  the  time 
step  size  (optimal  relaxation  factor)  to  minimize  the  time  steps  (number 
of  iterations)  for  obtaining  the  steady  state  converged  solution. 

The  NLMR  method  developed  by  Kennon  [7],  [8],  [9],  actually  can  be 
considered  as  a  time-dependent  method  with  variable  time  step  size.  The 
NLMR  method  is  based  on  the  sequential  minimization  of  the  global  resi¬ 
dual.  A  relaxation  factor  is  introduced  to  minimize  the  norm  of  the 
residual  at  each  iteration.  This  iteration-dependent  optimal  relaxation 
factor  drives  the  iterative  solutions  to  convergence.  In  order  to  obtain 
a  smooth  convergence  history,  several  options  are  available  to  this  ac¬ 
celeration  procedure.  For  instance,  each  three  regular  smoothing  iter¬ 
ations  are  followed  by  two  accelerated  iterations  [7].  Marchuk  [6]  has 
proved  that  under  certain  conditions,  the  variational  optimization  proc¬ 
ess  will  produce  the  highest  convergence  rate  for  iterative  method,  and 
the  norm  of  the  residual  will  form  a  monotone  decreasing  sequence  with 
Iteration.  Marchuk  [6]  also  pointed  out  that  the  convergence  speed  of 
the  minimum  residual  method  can  be  even  further  improved  by  using 
single-iteration,  two-step  minimum  residual  method. 

The  main  objective  of  the  present  study  is  to  extend  Marchuk's  idea 
to  generalize  the  NLMR  method  using  the  time-dependent  approach  and  the 
single-iteration,  multi-step  algorithm.  The  applications  of  the  GNLMR 
method  are  demonstrated  by  two  numerical  examples:  one  dimensional 
Burgers'  equation  and  the  two  dimensional  heat  conduction  equation.  Se¬ 
veral  interesting  problems  originating  from  this  method  such  as  inte¬ 
gration  by  sampling,  the  concept  of  grid-dependent  relaxation  factor,  and 
the  determination  of  the  global  minimum  are  also  discussed. 
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A. 3. 2  Theoretical  Aspects _ 


A.3.2.1  Optimization  of  the  Euler  Scheme  for  Linear  Problems 

Let  us  first  consider  a  well-posed  linear  initial  value  problem 


3$/3x  =  L$  -  F 

in 

2 

(1) 

*  =  *B 

on 

32 
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Applying  the  Euler  one-step,  time-consistent,  explicit  scheme,  the 
finite  difference  equation  of  (1)  can  be  written  as  : 

0t+1  =  "  f)  (4) 

where 

l  denotes  the  scheme-dependent  difference  analog  of  L, 
f  is  the  discrete  analog  of  F  and 
t  denotes  the  time  level. 

Most  linear  stability  analyses  of  the  scheme  represented  by  equation 
(4)  do  not  consider  the  effects  of  boundary  conditions,  thus  resulting 
in  overly  restrictive  and  even  incorrect  conclusions.  For  example,  H. 
D.  Thompson,  et  al .  [10],  proved  that  the  cell  Reynolds  number  re¬ 

striction  for  convection-diffusion  problems  derived  from  linear  stability 
analysis  is  a  commonly  accepted  misconception.  Moreover,  the  numerical 
experiments  performed  by  S.  R.  Kennon  [7],  [8],  [9],  using  the  NLMR  method 
showed  that  the  usual  Courant-Friedrichs-Lewy  (CFL)  number  limitation  for 
both  linear  and  nonlinear  problems  can  be  significantly  exceeded.  It 


should  be  pointed  out  that,  although  Thompson  et  al .  clarified  the  mis¬ 
conception  of  cell  Reynolds  number  limitation  for  1  inear 
convection-diffusion  problems,  their  results  did  not  give  the  best  value 
of  time  step  size  for  accelerating  the  scheme.  The  NLMR  method  provided 
a  simple  analytic  way  to  determine  the  optimal  acceleration  factors  for 
both  linear  and  nonlinear  problems.  However,  the  elementary  time  steps 
used  for  obtaining  the  corrections  still  follow  the  CFL  number  limitation 
concluded  from  the  linear  analysis. 

Assuming  that  coding  of  a  numerical  scheme  will  not  cause  too  much 
difficulty,  stability,  convergence  speed  (computer  time)  and  accuracy  of 
the  solution  are  three  major  factors  to  be  compromised  in  optimizing  the 
numerical  scheme.  We  now  present  an  easy  and  analytic  formulation  which 
allows  us  to  make  this  compromise  effectively. 

* 

Let  the  exact  solution  of  L$  -  F  =  0  be  denoted  as  $  . 

Define 


as  the  error  and  residual  vectors  at  time  level  t,  respectively.  Hence, 
the  time  evolution  of  both  error  and  residual  vectors  satisfy  the  fol¬ 
lowing  equations  : 


t  _ 

(7) 

,t+l 

=  rfc  ♦  Ax(ert) 

(8) 

.t+1 

» 

=  ?  *  AxCt^) 

(9) 
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The  residual  and  the  error  norms  at  time  step  t+1  can  be  expressed  as 


=  Ir'l’  ♦  ZAtfr*,  er1)  »  (Ax)1  liar'll*  (10) 

=  i?1)*  *  nl)  *  (ao’iu?1:1  (id 

Define  the  rate  of  convergence  T  and  rate  of  damping  A  by 

T  =  -log(  |lrt+*||/flrt||)  (3^) 

A  =  -log(!l5t+1il/!!Stil)  (13> 


The  convergence  and  stability  of  scheme  (4)  requires  that  both  T  and  A 
are  greater  than  zero.  However,  the  Lax  equivalence  theorem  [11]  states 
that  for  linear  initial  value  problems  ,  stability  is  the  sufficient  and 
necessary  condition  for  convergence.  Thus,  the  limitation  on  At  can  be 
obtained  by  solving  the  inequalities  : 

At  >  0  and  |!rt+1ii2  <  kV  (14) 
The  solution  of  the  above  inequalities  can  easily  be  obtained  as  : 

0  <  At  <  -2(r\  nrt)/;lertil2  (15) 

It  should  be  noted  that  equation  (15)  addresses  only  stability  and  not 
the  accuracy  of  the  solution.  Nevertheless,  if  time-accurate  solution 
is  required  (with  a  specified  accuracy),  equation  (15)  provides  the  free 
choice  of  At  to  meet  both  stability  and  the  desired  accuracy  conditions. 
On  the  other  hand,  if  time  evolution  is  not  important  ,  At  can  be  chosen 
as  large  as  possible  to  minimize  the  number  of  time  steps  for  obtaining 


the  steady  state  solution.  The  optimal  value  of  At  can  be  determined 
by  maximizing  the  convergence  rate  T  and  the  result  is 


(i0opt  =  -  (r\  (16) 

Thus,  for  steady  state  solution  of  equation  (1),  the  optimized  scheme 

*t+1  =  **  +  (At)opt(t^  -  f)  (17) 

will  produce  the  highest  rate  of  convergence.  The  monotone  convergence 
behavior  of  the  optimized  scheme  can  be  best  illustrated  pictorially  by 
figure  9  and  figure  10.  As  for  the  stopping  criteria,  the  L0  norm  or  L 

C.  • 

norm  of  the  residual  at  new  time  level  are  commonly  used  in  determining 

if  the  steady  state  solution  is  obtained.  However,  by  examining  equations 
(8)  and  (10),  it  is  interesting  to  note  that  At  will  approach  zero  as  the 
converged  solution  is  achieved.  This  provides  us  with  a  new  criteria  for 
stopping  the  iterative  process. 


A. 3. 2. 2  r  ti-Step  Minimum  Residual  Method  for  Linear  Problems 

The  method  described  in  the  previous  section  can  be  used  to  optimize 
the  single-iteration,  one  step  scheme.  However,  the  speed  of  convergence 
of  scheme  (17)  can  be  improved  even  further  by  multi-step  algorithm. 

Assume  that  M  steps  are  used  to  iterate  at  each  time  level.  Using  the 
Einstein  summation  convention  where  repeated  subscripts  are  summed,  the 
multi-step  algorithm  for  equation  (1)  is  then  defined  as  follows  : 
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+  u5  ;  m  =  1  ~  M 

m  m  * 


where 


6j  =  l<t>  -  f 


6m  =  *  ;  m  >  1 


are  the  residual  vectors  at  step  m,  and  u  are  the  relaxation  factors  to 

m 

be  determined  by  minimizing  the  l_2  norm  of  the  residual  at  time  level 
(t+1).  If  the  previous  definitions  of  error  vector  and  residual  vector 
are  used,  the  following  equations  can  easily  be  verified  : 


v  =  1C 

(20) 

t+1  =  rfc  +  u  26 
m  m 

(21) 

t+1  _  -t  .  . 

-  £  +  u  5 

mm 

(22) 

The  l_2  norm  of  the  residual  at  time  level  t+1  can  be  expressed  as 

«rt+1r  =  ||rV  +  2«m(r\  I6m)  +  (£6^  «„)»„,»„  :  m>  n  =  1  ~  M  (23> 


In  order  to  get  the  highest  rate  of  convergence,  clearly,  u>m  are  the 
solutions  of  the  following  system  of  linear  equations  : 


3l73w  =0  or  (rt,  25  )  +  (  25  ,  25)w  =  0 

m  m  m  n  n 


Multiplying  equation  (24)  by  wm  ,  it  follows  that 


Subtracting  equation  (25)  from  equation  (23)  and  using  equation  (24)  re¬ 
sults  in 


Thus,  the  residual  norms  for  the  multi-step  minimum  residual  method  also 
shows  a  monotone  convergence  behavior  which  guarantees  the  stability  of 
the  iterative  scheme  and  produces  the  highest  rate  of  its  convergence. 


A. 3. 2. 3  Optimization  of  the  Euler  Scheme  for  Nonlinear  Problems 

For  clarity,  we  consider  two-dimensional  problems  and  equations  in 
conservative  form  only.  The  extension  to  the  N-dimensional  problems  and 
non-conservative  equations  is  then  straightforward. 

The  conservative  form  of  the  governing  equations  for  most  engineering 
problems  can  be  written  as  : 

a*/3t  =  LyNv(0,  *x,  0y)  -  F  (27) 

where 

ll  =  3/3x  ,  L2  =  3/3y  (28) 

and  Nv  is  the  nonlinear  differential  operator  in  x  coordinates.  Usinq 

V  J 

the  Euler  one-step,  time-consistent,  explicit  scheme,  the  finite  differ¬ 
ence  form  of  equation  (27)  can  be  written  as  : 

0t+1  =  +  AxCe^U1,  ♦tx>  0ly)  -  f]  (29) 

or 
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r.on  -  -  -t -  - 


where 


0t+1  =  01  +  Air1 


r*  =  n^C^,  0tx,  0^y)  "  f 


is  defined  as  the  residual  at  time  level  t  . 

Therefore,  the  residual  at  time  level  t+1  can  be  expressed  as 


t+1  _  ,  .t+1  .t+1  .t+1  x  _  f 

r  =  *yN  (0  ,  0  x.  <P  y)  T 


After  expanding  the  nonlinear  operator  Nv  into  Taylor  series,  it  follows 


t^l  _  .  /  nV/  t  .  t  t  \ 

r  -  MN  (0  ,  0  x,  0  y) 


+  [(3Nv/0t)rt  +  (3Nv/30tx)(rt)x  +  (3Nv/30ty)(rt)y]  At  (33) 


,|V/o..t  w_t» 


+  0(At2)}  -  f 


In  summary, 


where 


*  rl  *  .„(AO- 


*j  =  tv[(3Nv/»t)rt  ♦  (3Nu/3^tx)(rt)x  ♦  (3Nv/3*ty)(rt)y]  C3S) 


and  the  coefficients  of  higher  order  terms  a^,  a^,  a^,  ...  can  be  de¬ 
termined  by  direct  Taylor  series  expansion.  Equation  (34)  indicates  that 
the  residual  at  time  level  t+1  is  a  polynomial  (hereinafter  called  resi¬ 
dual  polynomial  [7]  or  RP)  of  the  time  step  size.  The  total  number  of 
terms  in  the  residual  polynomial  (34)  depends  on  the  degree  of  nonline¬ 


arity  of  the  operator  N 
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If  N  is  polynomial  in  its  arguments  then  the  Taylor  series  truncates 
and  becomes  exact.  Thus,  the  l_2  norm  of  the  residual  at  time  level  t+1 

can  be  expressed  as 

||rt+1||2  =  ||rV  +  2(rt,  am)(AT)m  +  (am>  an)(Ax)m(AT)n  ;  m,  n  >  1  (36) 

Equation  (36)  implies  that  the  residual  norm  at  time  step  t+1  is  a  posi¬ 
tive  polynomial  (hereinafter  called  minimizing  polynomial  [7]  or  MP)  of 
the  time  step  size  (  to  be  determined  ).  Thus,  the  stability  of  scheme 
(29)  will  be  guaranteed  provided  that  At  is  chosen  as  the  optimizer  of 

the  minimizing  polynomial  (36)  such  that  flrt+^||  is  an  infimum  (global 
minimum).  However,  the  determination  of  the  optimizer  needs  special  nu¬ 
merical  techniques  [14].  The  rate  of  convergence  depends  on  the  relative 
difficulty  in  finding  the  optimizer.  To  reduce  this  difficulty,  the 

linearized  operator  of  Nv  may  be  applied.  If  Nw  is  truncated  to  the  first 
order  of  At  (linearized  operator),  the  approximate  residual  vector  is 

I 

(rt+1  )  =  r4  +  a2AT  (37) 

Then,  the  approximate  MP  is 

I 

l|rt+1  II2  =  IjrV  +  2( a x ,  rfc) At  +  (a^  a1)(AT)2  (38) 

The  optimizer  of  the  above  equation  can  be  easily  found.  It  should  be 
noted  that  equations  (36)  and  (38)  have  the  same,  but  negative  slope  at 
At  =  0  .  Figures  3  and  4  give  the  best  illustration  for  the  effects  of 
the  linearization  on  the  convergence  rate. 


A. 3. 2. 4  The  Generalized  Nonlinear  Minimum  Residual  Method 


The  GNLMR  method  actually  is  the  application  of  the  methods  described 
in  the  previous  sections.  Let  us  consider  the  problems  governed  by 
equation  (27)  and  assumme  that  M  steps  are  used  at  each  time  level  t. 
The  multi-step  algorithm  for  nonlinear  problems  is  defined  as 


,t+l 


+  um6m  +  °(u2m)  ;  m  =  1  ~  M 


(39) 


where  repeated  indices  are  summed. 
The  residual  at  step  m  is  defined  as 


5,  =  «vNV.  *\)  -  f 


(40) 


»„  =  Ij  ON',/3*t)61„_1  *  ON'’/S0tx)(6m.1)x  *  (5N''/a»ty)(6rn.1)y]  (41) 


The  coefficients  of  the  higher  order  terms  of  u  can  be  obtained  by 

m 


Taylor's  series  expansion.  If  only  linear  terms  of  wm  are  retained,  the 


residual  polynomial  (RP)  at  time  step  t+1  can  be  expressed  by  Taylor's 
series  expansion  as 


t+1  _  »  t+1  ,t+l  x  - 

r  =  *VN  (*  ,  <fi  x,  0  y)  “  f 


l  Nwi>t  +  «  6  ,  **  +  u  (6  )  ,  f  +  «  (6  )  ]  -  f 

v  1  m  m’  x  mv  irrx  y  nr  nryJ 


*  rt  *  tv([  ♦  (3NU/l»tx)(5m)x  ON'’/3,ty)(Jm)y  ]., 


(42) 
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Therefore,  the  minimizing  polynomial  (MP)  at  time  step  t+1  can  be  deter¬ 
mined  as 

l|rt+1 1|2  =  llrV  +  g(wm)  (43) 

where  g(«m)  Is  a  polynomial  in  um.  For  a  highly  nonlinear  differential 
equation,  g  will  be  a  complicated  multi-variable  polynomial  that  depends 
on  the  steps  we  used  and  the  degree  of  the  nonlinearity  of  the  differen¬ 
tial  operator  Nv.  Thus,  a  fast  and  accurate  method  to  determine  the  op¬ 
timizer  of  MP  is  required  for  the  GNLMR  method  to  guarantee  the  highest 

rate  of  convergence.  If  the  linearized  operator  of  Nv  is  applied  to  re¬ 
duce  the  difficulty  of  finding  the  global  minimum,  the  approximate  opti¬ 
mizer  of  (43)  can  be  determined  by  the  method  described  in  Section  2.2  . 

Since  the  coefficients  in  the  MP  are  obtained  by  integrating  the  res¬ 
iduals  over  the  whole  domain,  the  GNLMR  method  requires  a  large  amount 
of  computer  storage  to  save  the  residuals  from  each  step  m.  This  is  also 
the  common  problem  for  all  the  conjugate-gradient  type  methods  [12].  For 
new  generation  of  supercomputers  storage  may  not  cause  too  much  problem 
[12].  However,  if  storage  restriction  does  exist,  an  integration  by 
sampling  could  be  considered.  Nevertheless,  several  problems  accompany 
this  idea  [13].  Since  a  statistical  sampling  procedure  can  be  used  to 
get  the  random  data  for  integration,  the  question  of  error  incurred  by 
the  sampling  should  be  addressed.  If  the  error  caused  by  sampling  ap¬ 
proximation  can  be  determined,  it  should  be  possible  to  design  an  optimal 
sampling  procedure  to  minimize  the  error. 


The  problem  of  finding  the  global  minimum  (infimum)  or  the  global 
maximum  (supremum)  for  a  function  of  several  variables  is  a  basic  problem 
in  global  optimization  theory.  In  the  GNLMR  method,  the  optimal  relaxa¬ 
tion  factors  are  actually  the  optimizer  of  the  MP.  The  relative  diffi¬ 
culty  to  find  the  optimizer  depends  on  the  degree  of  non-linearity  of  the 
governing  differential  equations  and  the  steps  we  used  for  the  multi-step 
algorithm.  As  pointed  out  by  H.  Ratschek  and  J.  Rokne  [14],  interval 
method  is  the  only  among  all  the  existing  methods  that  guarantee  always 
the  location  of  the  global  extremum  with  arbitrary  accuracy.  Since  the 
interval  method  is  basically  iterative,  it  will  require  considerable 
computer  time  to  get  the  accurate  result  for  a  multi-variable  polynomial 
of  higher  order.  The  highest  rate  of  convergence  of  GNLMR  method  will 
be  guaranteed  if  the  optimizer  of  MP  can  be  efficiently  determined. 

Since  boundary  conditions  can  be  applied  exactly  in  the  regions  that 
are  neighboring  to  the  boundary,  it  is  obvious  that  during  the  process 
of  relaxation,  the  iterative  solutions  in  those  region  are  more  accurate 
than  the  regions  that  are  far  away  from  the  boundary.  This  implies  that 
the  relaxation  factor  should  vary  from  one  subregion  to  another,  or  even 
from  point  to  point  since  the  subregions  or  points  where  residuals  are 
large  definitely  need  more  correction. 

A. 3. 3  Numerical  Examples 

Two  numerical  test  cases  are  used  to  demonstrate  the  applications  of 
the  GNLMR  method  :  one-dimensional  Burgers's  equation  and  two-dimensional 
heat  conduction  equation.  Both  test  cases  were  computed  using 
non-accelerated  method  where  time  step  size  satisfies  the  CFL  limitation. 
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In  addition,  three  accelerated  computations  were  performed  with  M=l,  2, 
3  where  M  is  the  number  of  steps  used  in  the  multi-step  algorithm.  It 
should  be  pointed  out  that  when  applying  the  GNLMR  method,  it  is  not 
necessary  to  specify  the  time  step  to  obtain  the  residuals  at  each  step. 
However,  in  order  to  compare  the  relative  speed  of  convergence  between 
the  GNLMR  method  and  the  non-accelerated  method,  a  time  step  size  that 
satisfies  the  CFL  limitation  is  used  in  calculating  the  residual  at  each 
step  for  both  cases.  Moreover,  the  effect  of  linearization  to  the  rate 
of  convergence  is  investigated  by  solving  both  exact  and  approximate  MPs 
for  M=1  in  the  case  of  Burgers's  equation.  Comparisons  are  based  on  the 
accuracy  of  the  solution  and  the  actual  computer  time  required. 

A. 3. 3.1  Burgers'  Equation 

According  to  the  notations  defined  in  the  Section  2,  the 
one-dimensional,  viscous  Burgers'  equation  can  be  written  as 


where 


30/3t  =  LxN(0,  4>x) 

N(e,  0X)  =  -**/2  +  v*x 


(44) 


*1 


•.v. 


and  v  is  the  viscosity  coefficient.  In  this  example,  v  =  0.07  is  used. 
The  initial  and  the  boundary  conditions  are  chosen  as  follows 


♦(1,  O  =  0 

0(0,  t)  =  1 

*(x,  0)  =  1-x 


(45) 


The  FTCS  scheme  is  applied  to  discretize  equation  (44)  as 


♦  =  *  +  Ax{-[(^fci+1)1  -  (*  ,._1)i]/4Ax  + 

v(*  i  +  1  -  2/1  +  ^.jJ/Ax2} 


(46 


where  1  denotes  the  ith  grid  point  in  the  total  of  41. 
The  exact  RP  for  M=1  can  be  expressed  as 


where 


RP  =  t*  *  =  rt  +  a^Ai  +  a2(At)2 

al  =  a/axC-^r1  +  vr^) 
a2  =  -o.sa/axtr1)2 


(47] 


(48) 


The  optimal  value  of  At  is  chosen  as  the  optimizer  of  the  exact  MP 

MP  =  ||rt+V  =  AQ  +  A:At  +  A2(At)2  +  A3(At)2  +  A4(At)‘  (49) 

where 

A0  =  Hr1 1| 2  ,  Aj  =  2(r\  aj) 

a2  =  (ap  aj)  +  2(rt,  a2) 

A3  —  2(ap  a2)  ,  A^  -  (a2>  a2) 

If  linearized  operator  of  N  and  M  steps  are  used,  the  residual  polynomial 
is  truncated  up  to  its  first  order  as 


RP  =  rt+1  =  rl  +  a  u 
m  m 

yhere 

a_  =  a/ax[-^t5„  +  v(6jl 


(511 


and  6m  can  be  determined  from  equations  (40)-(41). 

The  minimizing  polynomial  MP  is  then 

MP  =  Ir1*1!’  =■  IrV  *  2(r\  »„)»„,  *  Um.  *.)■,»■ 


Thus,  the  optimal  relaxations  «m  can  be  easily  determined  by  solving  the 

following  system  of  linear  equations. 

Amnw  =  b 
mn  n  m 

where 


\in  ^am’  an^ 


b»  =  -<r  • 


and  A  is  a  symmetric  matrix  of  order  M. 
mn 

Figure  1  shows  the  exact  steady  state  solution  and  the  numerical  sol¬ 
utions  after  100  iterations.  It  is  obvious  that  the  GNLMR  method  gives 
the  most  accurate  results.  The  computer  time  costs  of  100  iterations  for 
each  case  are  shown  in  figure  2.  It  indicates  that  the  GNLMR  method  needs 
more  computer  time  per  iteration.  This  depends  on  the  number  of  steps 
that  were  used  in  the  multi-step  algorithm  and  whether  a  linearized  or 
exact  MP  was  solved.  However,  if  a  desired  accuracy  is  specified,  the 
GNLMR  method  needs  fewer  iterations  to  achieve  the  accuracy  requirement 
as  shown  in  figure  3.  Therefore,  the  overall  judgement  of  the  efficiency 
of  the  GNLMR  method  should  be  based  on  the  computer  time  rather  than 
number  of  iterations  required  to  obtain  a  solution  with  a  specified  ac¬ 
curacy.  Figure  4  does  prove  the  efficiency  gained  by  using  the  GNLMR 
method.  The  variations  of  the  relaxation  factors  with  respect  to  time 
(or  iteration)  show  that  a  sudden  change  in  their  magnitudes  occurred  just 
before  converged  solution  is  reached  as  shown  in  figures  5,  6  and  7.  It 
is  interesting  to  note  that  all  the  relaxation  factors  reduce  to  zero  as 
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converged  solution  is  obtained.  Since  the  relaxation  factor  is  equivalent 
to  the  time  step  size  in  our  formulations,  this  phenomena  can  be  inter¬ 
preted  as  that  the  time  marching  will  stop  once  the  steady  state  solution 
is  obtained.  This  provides  an  alternative  criteria  for  stopping  the  it¬ 
erative  process  in  the  GNLMR  method. 

It  must  be  mentioned  that  according  to  figures  3  and  4,  both  the  rate 
of  convergence  and  smoothness  of  convergence  of  the  NLMR  method  can  be 
improved  even  further  by  the  multi-step  algorithm.  Moreover,  figures  3 
and  4  show  a  peculiarity  that  during  the  first  few  tens  of  iterations 
solution  converges  very  slowly  and  then  suddenly  accelerate  to  its  maximum 
rate.  The  converged  solution  is  then  reached  almost  immediately  following 
the  maximum  rate  interval.  This  behavior  is  entirely  different  for  linear 
problems  as  mentioned  by  Marchuk  [6].  For  a  positive  definite  matrix  (a 
linear  operator),  the  convergence  rate  during  the  initial  iterations  is 
much  higher  than  the  asymptotic  rate  of  convergence  if  the  method  of 
minimum  residual  is  applied.  Hence,  in  practice  it  is  not  necessary  to 
solve  the  exact  MP  for  a  nonlinear  problem;  the  linearized  MP  can  be  used 
and  still  guarantee  a  high  rate  of  convergence  as  shown  in  figures  2,  3 
and  4. 

A. 3. 3. 2  The  Heat  Conduction  Equation 

The  unsteady  heat  conduction  equation  is  given  by 

3*/3x  =  Lo  ,  L  =  a(32/3x2  +  32/3y2)  (53) 

where  a  is  the  thermal  diffusivity. 
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and  A _  is  a  symmetric  matrix  of  order  M. 

mn 

The  numerical  results  are  summarized  by  figures  8  to  13.  As  indicated 
in  figure  8,  the  computer  time  required  per  iteration  for  the  GNLMR  method 
Is  longer  than  for  the  non-accelerated  method.  However,  figure  9  shows 
that  if  the  maximum  number  of  iterations  is  limited  to  80,  the  accuracy 
that  can  be  obtained  by  the  GNLMR  method  is  better  than  that  of  the 
non-accelerated  method.  If  a  fixed  computer  time  is  specified,  figure 
10  clearly  shows  that  the  GNLMR  method  will  produce  the  best  solution. 
Figures  11  to  13  show  the  variations  of  the  optimal  relaxation  factors 
versus  number  of  iterations.  Figures  9  and  10  confirm  that  using  the 
multi-step  algorithm  improves  both  the  smoothness  and  the  rate  of  con¬ 
vergence  for  the  NLMR  method.  If  the  GNLMR  method  is  applied  to  a  linear 
problem,  it  is  interesting  to  note  that  the  rate  of  convergence  during 
the  initial  few  iterations  is  much  higher  than  the  asymptotic  rate  of 
convergence  as  indicated  by  figures  9  and  10. 

It  should  be  pointed  out  that  the  linear  version  of  the  GNLMR  method 
is  similar  to  the  generalized  minimum  residual  GMRES  method  developed  by 
Saad  and  Schultz  [15]  to  solve  non-symmetric  linear  systems  of  equations. 
This  method  was  recently  modified  by  Wigton  et.  al.  [12]  to  solve 
non-linear  problems  in  gas  dynamics  using  Newton's  iteration  method. 
Since  Gram-Schmidt  orthogonal ization  procedures  are  applied  to 
orthogonal i ze  the  search  directions  at  each  iteration  in  the  GMRES  method, 
this  method  demands  a  large  number  of  arithmetic  operations  and  a  large 
computer  memory.  Moreover,  as  pointed  out  by  Marchuk  [6],  the  implemen¬ 
tation  of  the  orthogonal i zation  algorithm  with  respect  to  high-order  ma¬ 
trices  usually  results  after  a  few  tens  of  iterations  in  a  numerical 
Instability  because  of  nonlinearity.  This  implies  that  the  maximum  number 


of  steps  which  can  be  used  in  one  iteration  for  the  multi-step  algorithm 
is  practically  limited  to  a  moderate  finite  number. 

A. 3. 4  Conclusion  Remarks 

The  GNLMR  method  and  its  applications  were  presented.  The  accelerating 
mechanism  and  the  monotone  convergence  behavior  of  the  GNLMR  method  has 
been  proved  by  theoretical  studies  and  the  numerical  experiments.  It  was 
found  that  both  the  rate  and  the  smoothness  of  convergence  of  the  NLMR 
method  can  be  improved  even  further  by  the  optimized  multi-step  algorithm. 
Both  numerical  test  cases  confirmed  that  all  the  optimal  relaxation  fac¬ 
tors  vanish  as  converged  solution  is  achieved  thus  providing  a  new  stop¬ 
ping  criteria  for  the  iterative  process.  The  numerical  results  also 
proved  that  linearization  of  a  nonlinear  operator  in  the  GNLMR  method 
reduces  the  problems  of  determining  the  optimizer  of  the  MP.  This 
linearization  still  guarantees  a  high  rate  of  convergence.  Hence,  the 
GNLMR  method  can  be  simplified  in  practical  problems  when  iteratively 
solving  nonlinear  partial  differential  equations. 
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Equations",  AIAA  Paper  86-0554,  January  1986. 

4.  Huang,  C.  Y.,  Kennon,  S.  R.  and  Dulikravich,  G.  S.,  "Generalized 
Non-Linear  Minimal  Residual  (GNLMR)  Method  for  Iterative  Algorithms", 
submitted  to  the  Journal  of  Computational  and  Applied  Mathematics. 

A. 5  List  of  Professional  Personnel 

Dr.  George  S.  Dulikravich  —  Assistant  Professor,  ASE/EM  Department 

Stephen  R.  Kennon  —  Ph.D  Candidate  In  ASE/EM  Department 

Chung-Yuan  Huang  —  Ph.D  Candidate  in  ASE/EM  Department 

A. 6  Interactions  (Coupling  Activities) 

Papers  presented  at  technical  meetings  : 

1.  "Optimal  Acceleration  Factors  for  Iterative  Solution  of  Linear  and 
Non-Linear  Differential  Systems"  by  S.  R.  Kennon  has  been  presented 
at  the  AIAA  23rd  Aerospace  Sciences  Meeting,  Reno,  Nevada,  January 
14-17,  1985. 
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2.  "Direct  Solution  of  Polynomial  Systems  of  Equations  With  Applications 
to  Euler  Equations",  will  be  presented  by  S.  R.  Kennon  at  AIAA  24th 
Aerospace  Sciences  Meeting,  Reno,  Nevada,  January  6-9,  1986. 


Consultative  and  Advisory  Functions  : 

Dr.  G.  S.  Dulikravich  delivered  the  following  invited  lectures  that  in¬ 
cluded  the  optimal  acceleration  factor  theory  : 

—  Mech.  Eng.  Dept.,  U.  of  California,  Davis,  CA,  Jan.  1985. 

—  Mech.  Eng.  Dept.,  Rice  Univ.,  Houston,  TX,  Jan.  1985. 

—  General  Electric  Co.,  Evandale,  OH,  April  1985. 

—  Aerospace  Eng.  Dept.,  U.  of  Texas,  Arlington,  TX,  May  1985. 

—  Mech.  Eng.  Dept.,  Old  Dominion  Univ.,  Norfolk,  VA,  May  1985 
--  Allison  Gas  Turbines,  Indianapolis,  IN,  July  1985. 

—  Mech.  Eng.  Dept.,  U.  of  Texas,  Austin,  TX,  Oct.  1985. 

—  Aerospace  Eng.  Dept.,  Pennsylvania  State  U.,  Univ.  Park,  PA,  Nov. 
1985. 

—  Mech.  Eng.  and  Material  Sciences  Dept.,  Duke  U.,  Durham,  NC,  Dec. 
1985. 

Mr.  S.  R.  Kennon  delivered  the  following  invited  lectures  that  included 
the  optimal  acceleration  factor  theory  : 

—  Dept,  of  Mathematics,  U.  of  Texas,  Austin,  TX,  March  1985. 

--  Computational  Fluid  Dynamics  Branch,  NASA  Ames  Research  Center, 
Field,  CA,  August  1985. 

—  Computational  Aerodynamics  Branch,  Lockheed  -  California  Co., 
Burbank,  CA,  August  1985. 
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A  mini-grant  was  obtained  from  the  Computational  Fluid  Dynamics 
Branch  at  NASA  Ames  Research  Center  (  $7,000  +  5  hours  of  CRAY  X-MP  com¬ 
puter  time  )  that  gave  an  opportunity  to  Mr.  S.  R.  Kennon  to  visit  there 
last  summer  for  six  weeks  and  to  develop  the  basic  algorithm  for  his 
latest  AIAA  paper  on  polynomial  systems  of  equations  for  fast  iterative 
algorithms. 

A. 7  New  Discoveries,  Inventions  or  Patent  Disclosures 
No  contribution. 

A. 8  Additional  Information  for  Evaluation  of  the  Progress 

Work  in  progress  is  oriented  towards  accelerating  implicit  iterative 
methods  for  the  Euler  and  Navier-Stokes  equations  of  gasdynamics.  Future 
work  will  include  (a)  accelerating  the  solution  of  the  Euler  and 
Navier-Stokes  equations  and  (b)  further  exploitation  of  the  polynomial 
properties  of  the  Euler  and  Navier-Stokes  equations  by  applying  classical 
algebraic  polynomial  elimination  theory. 

The  immediate  task  is  an  application  of  the  generalized  non-linear 
minimal  residual  method  to  the  finite-volume  Runge-Kutta  time-stepping 
scheme  of  Jameson.  This  scheme  uses  multiple  time  steps  that  are  pres¬ 
ently  limited  by  the  linear  stability  analysis.  It  is  expected  that  our 
method  of  multiple  optimal  acceleration  factors  can  accelerate  this  pop¬ 
ular  scheme  and  enhance  its  stability. 
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B.  Report  of  the  Group  Headed  by  Or.  D.M.  Young,  Jr. 

B .  1  Abstract 

A  composite  adaptive  procedure,  based  on  the  use  of  variational 
techniques,  has  been  developed  for  the  automatic  determination  of  split¬ 
ting  parameters  for  accelerated  iterative  procedures.  Examples  include 
the  symmetric  SOR  method  and  the  shifted  incomplete  Cholesky  method,  each 
with  conjugate  gradient  acceleration. 

B.2  Research  Objectives  and  the  Statement  of  Work 

Mr.  Mai,  working  with  Professor  Young,  has  been  developing  and 
testing  "composite"  adaptive  procedures  for  speeding  up  the  convergence 
of  certain  iterative  algorithms  for  solving  large  sparse  linear  systems 
of  the  type  arising  in  the  numerical  solution  of  partial  differential 
equations.  The  iterative  algorithms  considered  include  a  basic  iterative 
method  combined  with  the  acceleration  procedure,  such  as  Chebyshev  ac¬ 
celeration  or  conjugate  gradient  acceleration.  The  effective  application 
of  these  algorithms  may  require  the  choice  of  "splitting  parameters"  as¬ 
sociated  with  the  acceleration  procedure. 

In  the  past  these  parameters  have  been  determined  adaptively,  but 
the  determination  of  the  splitting  parameters  has  been  done  separately 
from  the  determination  of  the  acceleration  parameters.  A  "composite" 
adaptive  procedure  is  designed  to  determine  both  sets  of  parameters  si¬ 
multaneously.  The  procedures  used  are  based  on  the  application  of  certain 
variational  principles  combined  with  the  solution  of  certain  nonlinear 
optimization  problems.  Preliminary  testing  has  been  done  on  generalized 
conjugate  gradient  acceleration  procedures  applied  to  the  symmetric  sue- 
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r. 
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B.6  Interactions  (Coupling  Activities) 

It  is  planned  to  describe  the  work,  in  a  paper  to  be  presented  at  the 
International  Conference  on  Vector  and  Parallel  Computing,  Leon,  Norway, 
June  2-6,  1986. 

B.7  New  Discoveries,  Inventions  and  Patent  Disclosures 
No  contribution 


B.8  Additional  Information  for  Evaluation  of  the  Progress 
No  contribution 
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